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o. 

The ergodic hypothesis asserts that a classical mechanical system will in time visit every available 
configuration in phase space. Thus, for an ergodic system, an ensemble average of a thermodynamic 
quantity can equally well be calculated by a time average over a sufficiently long period of dynamical 
evolution. In this paper we describe in detail how to calculate the temperature and chemical potential 
from the dynamics of a microcanonical classical field, using the particular example of the classical 
modes of a Bose-condensed gas. The accurate determination of these thermodynamics quantities is 
essential in measuring the shift of the critical temperature of a Bose gas due to non-perturbative 
many-body effects. 
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I. INTRODUCTION 



The shift in critical temperature of the homogeneous Bose gas has been the subject of numerous investigations 
over the past fifty years. As the density of this idealized system is constant, the shift due to the mean-field is zero, 
and the first order shift is due to long-wavelength critical fluctuations. The first estimates were due to Lee and Yang 
P, Q, who gave two different results for the dependence on the s-wave scattering length a. In 1999 Baym et al. 3] 
Q ' determined that the result should be 

AT C /T C ° - can 1 / 3 , (1) 

where n is the particle number density, and c is a constant of order unity. Several authors have attempted to calculate 
■^|- ■ this constant, and a wide range of results have been obtained, as summarised in Fig. 1 of [4]. However recent Monte 
£**>«. ' Carlo calculations have apparently settled this matter, giving a combined result of c m 1.31 ± 0.02 Hi]. A useful 
^sO , summary and discussion of this topic is provided by Andersen @ and Holzmann et al. 

Previously we have performed numerical simulations of an equation known as the Projected Gross-Pitaevskii equa- 
tion (PGPE), which can be used to represent the highly occupied modes of Bose condensed gases at finite temperature 
[H, S floL fill . [HJ . This equation is observed to evolve randomised initial conditions to equilibrium for which it is pos- 
sible to measure a temperature [l3|, [IHdH. The PGPE is non-perturbative, and hence includes the effect of critical 
fluctuations. The only approximation made is that the modes of the gas are sufficiently highly occupied as to be 
well described by a classical rather than a quantum field. The occupation condition is that mode k must satisfy 
I i (JVfc) ^> 1; however for practical simulations we may choose, for example, (Nk) > 5 [16]. This method is suitable for 
investigating many problems of current interest in ultra-cold Bose gases, such as the shift in critical temperature due 
Ch , to non-perturbative effects [17| . 

The PGPE describes a microcanonical system, with the classical field restricted to a finite number of modes for 
which the occupation number condition is met. In order to study the problem of the shift in T c it is necessary to 
accurately determine thermodynamic quantities defined as derivatives of the entropy such as the temperature and 
chemical potential. In 1997 Rugh developed a new approach to statistical mechanics where he derived an expression 
from the Hamiltonian of a system for which the ensemble average gives the temperature within the microcanonical 
ensemble [181 ] . However, if the system is known to be ergodic then the equilibrium temperature can be determined 
from the system dynamics over a sufficiently long period of time. 

We have applied an extension of Rugh's method to the PGPE Hamiltonian, and the appropriate expression to 
determine the temperature is given as Eq. (22) of [Hj]. This method was found to agree with the less rigorous 
methods described in [l3|, [l4j • In [l!| we made use of this method to calculate the shift in the critical temperature of 
the homogeneous Bose gas. Despite the calculation being performed with limited statistics and suffering from finite 
size effects, it gave a result of c = 1.3 ± 0.4 in agreement with the Monte Carlo results [1, Q. In [TtJ we applied this 
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method to the experiment of Gerbier et al. [20] who measured the shift in critical temperature of a trapped Bosc 
gas, and found good agreement with experiment. In this paper we give the details of our implementation of Rugh's 
method for a general mode basis for the PGPE. 



II. BACKGROUND 
A. Formalism of Rugh 

We consider a classical system with M independent modes. The Hamiltonian can be written as H = H(T), where 
r = jr.,} = {Qi, Pi} is a vector of length 2M consisting of the canonical position and momentum co-ordinates. We 
define the gradient operator V in terms of its components V, = d/dTi. In the notation of Rugh [2l[, the Hamiltonian 
H may have a number of independent first integrals, labelled by F = F\ , . . . , F m , that are invariant under the dynamics 
of H . A particular macro-state of such a system can be specified by the values of the conserved quantities, labelled 
as H = E,Fi = h. 

The usual expression for the temperature of a system in the microcanonical ensemble is given by 



1 fdS 



k B T \dE / 

where all other constants of motion are held fixed, and where the entropy can be written 



(2) 



e 



= jdT 6[E H(T)] HSih Fi(T)]. (3) 



Using Rugh's methods, the temperature of the system can be written as 

1 



P-Xr(r) , (4) 



k B T 

where the angle brackets correspond to an ensemble average, and the components of the vector operator T> are 

= (5) 

where can be chosen to be any scalar value, including zero. The vector field Xj- can also be chosen freely within 
the constraints 

VH ■ X T = 1, T>Fi ■ Xy = 0, 1 < i < m. (6) 

Geometrically this means that the vector field has a non-zero component transverse to the H = E energy surface, 
and is parallel to the surfaces Fi — li . The expectation value in Eq. ^ is over all possible states in the microcanonical 
ensemble; however if the ergodic theorem is applicable then it can equally well be interpreted as a time-average. For 
further details on the origin of this expression we refer the reader to Rugh's original papers [l8l l2ll . [22| . as well as 
derivations found in Giardina and Levi [2^|, Jepps et al. [24[ and Rickayzen and Powles [25| . 



B. Dimensionless Hamiltonian 



The classical Hamiltonian for the dilute Bose gas in dimensionless form is 

Cnl 



H 



^*(x)fr„ P v>(x) + v ad (x)|v(x)|' 



l^(x)| 4 



(7) 



where H = H / (Nc^l), Nc is the number of particles in the system, x = x/L, L is the unit of length, = h 2 / (2mL 2 ) 
is the unit of energy, and m is the mass of the particles. The dimensionless classical Bose field operator V'(x) is here 
normalized to one, J d 3 x|i/>(x)| 2 = 1, and C n \ is the nonlinear constant defined as 



Cnl 



N C U 8iraN c 



L 3 



(8) 
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where Uq — A-Kh 2 a/m. H sp is the single particle Hamiltonian with eigenbasis if sp 0fc(x) = £fc0fc(x), and I4d(x) is the 
dimensionless version of any potential additional to that included in the single-particle Hamiltonian. 

We restrict our field ^>(x) to be formed from the modes of the classical region C defined by a high energy cutoff in 
the single-particle basis, such that we can write 

V(x) = J2 <*4>k(*)- (9) 
The equation of motion for this restricted system is the Projected Gross-Pitaevskii equation 

,C^(X) = H ^ + p{ y (x) ^ (x)+C - |^ (x) |2^ x)}- (10) 
OT 

The projection operator V{F} projects function F onto the classical modes of the single particle Hamiltonian H sp 
via 

V{F( X )} = ]T 0*(x) / d 3 x' ^(x')^(x'). (11) 

It eliminates collisional processes which would cause the scattering of particles into higher energy single particle modes 
not represented by the classical field. 

In the single-particle basis of H sp the Hamiltonian (J7J) can be written 

H=J2 e nC* l c n + y mn c* m c n + -^- ^ c* m c* n c p c q (mn\pq) , (12) 

where we have defined the matrix elements 

V mn = J d 3 xft(x)y ad (x)^(x), (13) 
(mn\pq) = J d 3 x^(x)C(x)</) p (x)0 9 (x). (14) 
We make use of the real, canonically-conjugate coordinates Q n and P„ defined by 

« + c„), P n =n/%«-c ri ), (15) 



V2^ V 2 

with the inverse transformation 



€ f Qn + ^ =Pn) C * n = J'fQn- ~i= Pn. (16) 
Z V ^Cn V ^ V -^n 



C. Choice of vector field X 



In order to satisfy the conditions © we can choose a vector field of the form 

m 

X T = aVH + YJ& i X>F i , (17) 

i=l 

where the m + 1 coefficients {a,&i} are determined by the m + 1 simultaneous equations in Eq. ©. Due to the 
freedom in the choice of the vector operator T> we can set any number of components of the length 2M vector Xy 
to zero. This turns out to be useful as the components corresponding to the momentum and position variables can 
be different orders of magnitude. Two particular choices we make use of later are Xx,p with V = Up = {0, d/dPi} 
and Xt,q with V = Vq = {d/dQi, 0}. These lead to two different calculations for the temperature that only agree 
in general if the system is in thermal equilibrium. This provides a useful check that the numerical simulations have 
indeed reached thermal equilibrium, as well as giving two independent values for the temperature. 

For the classical Bose gas Hamiltonian ([7]) there can be several constants of motion that must be taken into account. 
These depend on the details of the single-particle Hamiltonian H sp — examples include components of the angular 
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and/or linear momentum. The effect of including these additional first integrals in the definition of the vector field 
is to account for the energy that is associated with a conserved quantity and hence is unavailable for thermalization. 
This ensures that only the appropriate free energy is used to calculate the temperature. We conjecture that the same 
result can be achieved by first transforming to a co-ordinate system where the total angular and linear momenta, etc, 
are all zero and therefore do not contribute to the energy of the system. Rugh demonstrated this explicitly in [21j for 
a system of particles with a conserved centre-of-mass motion, and our numerical results support this conjecture. 

An exception to this, however, is the conservation of normalization N = c*c„. This must be considered explicitly 
because there is no co-ordinate system in which it can be made to vanish. The constraint on TV means that the ground 
state of the system will, in general, have a finite energy that is not accessible for thermalization. We note that the 
effect of this constraint is in general more complicated than a simple subtraction of the ground state energy (which 
could be achieved by hand) and depends on the definition of the operator V used to calculate the temperature, as 
shown below. 

We need to choose a vector field Xj- which satisfies Eqs. ([6]) with F\ = N = c n c n- The result is 

x T = ™-x„vn (18 ) 
\T>H\ 2 - \ N VN-T>H' K ' 



where the parameter 



VN -VH 

Xn \vn\ 2 ' (19) 



looks similar to a chemical potential. Substituting Eq. (fTg)) into Eq. 0$ we find that our full expression for the 
dimensionless temperature T = hsT/Nc^L is 



1 _ / V 2 H - \ N V 2 N - V\ N -VN\ I {VH - \ N VN) ■ [V\VH\ 2 - (VH ■ VN)V\ N - \ N V{VH ■ VN)} \ 



T 
where 



VH\ 2 -\ N (VH -VN) / \ [\VH\ 2 - \ N (VH -VN)] 2 / 



,(20) 



VX N = ' + —^^-V(\VN\ 2 ). (21) 



V(VN ■ VH) VN ■ VH 
\VN\ 2 + \VN 

The chemical potential of the system can be determined in a similar manner using a different choice of vector field 
In particular we wish to calculate 



M _ ( dS\ 
k B T \dN 



E 



^•x M (r)), (22) 



where the conditions on the vector field are 

VH ■ X M = 0, VN ■ X p = 1. (23) 

This results in expressions identical to the RHS of Eqs. (|18ll9l20p but with H and N interchanged, and so we will 
not write them out in full. 

The technique outlined above has been described previously by one of us in Ref . , and in particular the result of 
Eq. (|20p was given in that work. The numerical method has also been applied by the current authors in Refs. [3, [13] • 
The purpose of this paper is to outline in detail the procedure for calculating the terms in Eqs. (|20|) and ([21) . This 
is not a trivial exercise, and the proliferation of terms can make it difficult to avoid errors in both the analytical 
and numerical implementation. By discussing our approach to evaluating Eqs. (|20]) and (|2"Tj) we hope to facilitate 
other researchers wanting to make use of the method. The basic numerical requirement is an efficient and accurate 
numerical transform from real space to basis space and vice versa. This may be provided, for example, by fast Fourier 
transforms using a basis of plane waves appropriate for a homogeneous system jUGS, 

or an efficient quadrature for 

trapped systems [l5j . 



III. FORMULAE 



In this section we give explicit analytic expressions for all quantities required in the evaluation of Eq. (|20|) and 
PTjl . As we are using the canonical co-ordinates defined in Eq. ([15]) it is easiest to evaluate the derivatives of the 
Hamiltonian using the mode expansion of the Hamiltonian (|12l) . However, this leaves us with expressions involving 
inner products of vector quantities with summations that would be computationally expensive to evaluate directly 
for any sizeable basis. We show how to define auxiliary field functions that can be used to simplify these terms for 
implementation with an efficient numerical transform. 
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A. Derivatives 



We begin by writing out the necessary derivatives of the Hamiltonian using the two choices of vector operator 
defined earlier, T>q and T>p. In component form the first derivatives of the Hamiltonian are 



{T> Q H)i = e? Q l + J j i v m c n + VtnCn] + C n u ^ [c*c*c m (pq\mi) + c* m c p c q {mi\pq)] 



(VpH)i = P t + —== ^ [Vmc* - V in c n ] + i—j== Y [ c p c *q c ™(Pl\ mi ) ~ c* n c p c q (mi 



2c, 



and the second derivatives are 



(T>QH)ij = Sijtj + —^-[V^ + Vji) + -^-^/elej ^ [ic* p c q (pi\qj) + c p c q (ij\pq) + c*c*(pq\ij}] 



(V 2 P H) t , = % + 5-=^ + V^} + 



2 
Cni 



: [4c*c q (pi\qj) - c p c q (ij\pq) - c*c q (pq\ij) 



(24) 
(25) 

(26) 
(27) 



We also need the derivatives of the other first integral of the Hamiltonian — the the normalisation given by N 
Y,kec \°k\ 2 - These are 



{V Q N)i = e l Q l , 
(V 2 Q N) ij =e i 5 ij , 

From Eq. (|20|) we can see that we need the terms 
V{\VH\ 2 )i = 2^(X>iJ)j(I? 2 iJ) 



(VpN)i = 
(DpN)ij = § Ji, 



V(\VN\ 2 ) l = 2Y J (DN) J {V 2 N) lJ1 



(28) 
(29) 

(30) 



V{VN ■ VH)i = [i'DN) ] (V 2 H) t3 + (Pff), (P 2 ^)* 

3 

Some of these expressions are simple to calculate e.g. 

V Q (\V Q N\% - 2]T(P Q iV) J (2?|A%- 

3 

= 2j> J Q i )(e i * ii ) = 2e 2 Q t , 



(31) 

(32) 
(33) 



VpdVpN] 2 ), = 2Y,(ppN) j {V P N) i 



(34) 
(35) 



(36) 



(37) 
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However, the remainder are more complicated. To calculate them efficiently we make the following vector definitions 

(PpH)j 



a 3 = V^( V Q H )h b 3 

fi = Vej(1>QN)j, 9 3 
with the corresponding auxiliary field functions 



A ( x ) = £oj0j( x ). b ( x ) = £Mj( x ). 
F(x) = 2 /^,(x), G(x) = ^ ^(x). 



(38) 
(39) 

(40) 
(41) 



Making use of these definitions we find 
V Q (\V Q H\% = 2^(P Qj ff) j (©2 F) . 



= *E 



1 " ••.<''.' " — + Vji] + ^V^7 E [ 4c p c «(^kj) + c p c q (ij\pq) + c*c*(pq\ij)] ) , 



= 2aief i^^aj^j + V^] + Cni^E^ [4c£c 9 (2%j) + c p c q (ij\pq) + c* p c* q (pq\ij)] , (42) 
V P (\V P H\ 2 ) t = 2Y,{VpH) j (V 2 P H) ij , 

3 

= 2 J2V^3 b 3 I % + 7= Wi3 + ^jfi] + 9-1= E [ 4c p c <jH4?> - CpC q {ij\pq) - c* p c* q {pq\ij)} j , 

1 C 

= 2^ + -= 2 6^ + V^-i] + ^ &j [4tjc,(pi|«) - c p c q {ij\pq) - c* p c* q (pq\ij)] , (43) 



and 



3 3 v \ p<j , 



= /*c 2 + y E + ^ + -f v^E [h^h^ + + c ; c *<^i^->] , (44) 

3 ' pqj 

Y,V>P*0iV>pH)ii = E V^9j (iij + 7T^=[Vi3 + Vji] + E [^* p c q {pi\qj) - <h>c q (ij\pq) - c* p c* q {pq\ij)] j , 

j j \ V J V J pq / 

= V^Ift + — = E 9j[Vij + Vji] + —y= E 9j [^c* p c q {pi\qj) - c p c q {ij\pq) - c* p c* q {pq\ij)] . (45) 



B. Evaluation of terms with auxiliary fields 

The above expressions are written in a form that they can be easily calculated using efficient quadratures in real 
space. We outline in detail how to calculate all of these terms, and write them in the most efficient form we have 
found. We make the definitions 

E^« c « = /d 3 x^*( x )Kd(x)^(x) = W u (46) 
E y ™< = / d 3 x^*(x)T4 d (x)0 4 (x) = W* , (47) 



where we have assumed that the potential V a( j(x) is real. We can therefore write the appropriate parts of Eqs. (|24|) 
and (l2"5l) as 



The next set of terms are 



X)[Ki<+^„c„] = 2Re(WU (48) 

n 

Y,[ v m<-V in c n ] = -2iIm(Wi). (49) 
Y,c* m c p c q (mi\pq) = /'d 3 x0*(x)^(x)|V(x) = K h (50) 

pqm 

Y,c*Ac m (j>q\mi) = /d 3 x|^(x)|V(x)</»i(x) = K*. (51) 

Hence, parts of Eqs. ([24]) and ([25]) can be written 

X] [ c t c *q c m(pq\mi) + c* n c p c q (mi\pq)] = 2Re(^), (52) 

pqm 

^ [c*c*c m (pq\mi) - c* m c p c q {mi\pq)] = -2ilm{Ki). (53) 

pqm 

For the terms involving the auxiliary fields we have expressions like 

X> iC ; C? (f%i) = / d 3 x^(x)|^(x)| 2 yl(x) = (AxJj, (54) 
pqj 

^2a jCp c q {ij\pq) = J d 3 x0*(x)A*(x)7/;(x) 2 = [A 2 ) h (55) 

^2a jC ;c* q (pq\ij) = J d 3 x^*(x) 2 <Mx)A(x) = (A 2 )*, (56) 

$>,K, = J d 3 x0*(x)Kd(x)A(x) = (A 3 ) i5 (57) 
j 

J^ajVji = J d 3 xA*(x)T4d(x)^(x) = {A s )*, (58) 
i 

where we have made use of the fact that Oj — a* . We have similar expression for all the other auxiliary vector fields 
b, /, and g and will refrain from writing these all out. We can write parts of Eqs. (|4"2"]) . (|4"3"]l . (|4"4"]) and (j4"S"|) as 

[4c;c,(pi|9j)+cpc,(ij|p9)+c; c ;(pg|ij)] = 4(Ai)< + 2Re(A2)i, (59) 

bj [4c; Cq (pi\qj) - c p c q (ij\pq) - c* pC ;(pq\ij)] = 4( J B 1 ) i - 2Re(5 2 ) l , (60) 

Pqj 

]T /, [4c;c q (pi\qj) + c p c q (ij\pq) + c* p c;(pq\ij)] = 4(*i)< + 2Re(F 2 ),, (61) 
pqj 

Y^9j [^c q {pi\qj) -CpCgiijlpq) - c p c*(pq\ij)] = 4(Gi)( - 2Re(G 2 ) i! (62) 

X a, [VJj + Vji] = 2Re{A 3 ) ll £ 6, + ty] = 2ReOB 3 )i, (63) 
E^[^ + ^] = 2 MF 3 h J^dAVij+Vji] = 2Re(G 3 ) t . (64) 
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The remaining terms are diagonal second derivatives of the form T> 2 H = '}2 i {T>' 2 'H)u, in which we have matrix elements 

J2c;c q (pi\qi) = /d 3 x|^(x)| 2 ^(x)| 2 = (Mi)*, (65) 

pq J 

^ Cp c 9 (u| M ) = /d 3 x0*(x)V(x) 2 = (Ma)*, (66) 

P9 J 

Y,^*(j>q\ii) = j d?^* = (M 2 )*, (67) 

and thus the appropriate parts of Eqs. ([26]) and |27|) are 

£ [Ac* p c q (pi\qi} + c p c q (ii\ P q) + c; c * q (pq\ii)] = 4(M 1 ) i + 2Re(M 2 ) i , (68) 

pq 

^[Aclc q {pi\qi) - c p c q (ii\pq) ^ c* p c* q {pq\ii)] = 4(Mi)j - 2Re(M 2 ) l . (69) 

We have now explicitly outlined the calculation of all terms necessary to calculate the temperature and chemical 
potential of the equilibrium dynamical evolution of a microcanonical PGPE system. 

C. Summary 

Here we summarise all the components of the vector quantities required for the evaluation of the temperature 
according to Eq. (|20|) making use of the definitions made in the previous section for both choices of vector operator 
T>q and Dp. The inner products between these vectors found in Eq. (|20|) can be easily evaluated numerically. In 
addition, these are all the terms required for calculating the chemical potential as well. 

{VqH), = e 2 Q l + V2^Re{Wi) + C„i V / 2^Re(^), (70) 

(P P H)i = Pi + ^-lmtWi) +C n Hpm(if i ), (71) 

V e; 

(V 2 Q H) U = e 2 + e 4 ^ + C nl e l [2(M 1 ) l +Re(M 2 ). i ], (72) 

(D P H) U = l + Y^ + ^[2(M 1 ) i -Re(M 2 )i}. (73) 

VqUVqNW = 2e 2 Q 4 , (74) 
Vp(\V P N\% = 2§, (75) 

J2CDQH)^QNh = e % {V Q H) l7 (76) 
j 

3 

V Q {\V Q H\% = 2a l e 3 l /2 + 2^r i Re(A 3 ) t + C ni ^-mA 1 ) l + 2Re(A 2 ) l }, (78) 
Vp{\V P U\% = 2 v ^ + -^Re(B 3 )< + -£[4(Bi)i-2Rfi(B 2 )<], (79) 



X^qJVJ^JI)^ - / l e- /2 + V7^Re(F 3 ) l + C nl V^[2( J F 1 i) t +Re( J F 1 2 ) 4 ], (80) 

i 

Y / (^pN) J (V P H) i] = ^ + -LHe(G 3 )i + ^[2(Gi)i - Re(G 2 ),]. (81) 
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IV. NUMERICAL RESULTS 

In this section we outline the algorithm for applying the formulae of the previous section to simulation results, and 
give examples of the results obtained for the PGPE for trapped Bose gas systems. 

We have previously demonstrated in [l3|, [TJ, [l5| that the PGPE will evolve any randomised initial field function 
to an equilibrium that can be characterised by a temperature, and have made use of this in the beyond mean-field 
calculation of the critical temperature of a trapped Bose gas in [l7| • 

The basic requirement to determine the temperature in such calculations is to evaluate Eq. (|20p over a sufficiently 
long period of simulation time for the ergodic averaging to be effective. Our procedure is as follows. Beginning with 
an initial random field configuration with a predetermined energy, we run our simulation for sufficiently long such 
that we can be reasonably sure that equilibrium has been reached (this can be checked as described below). This is 
typically for one hundred or more trap periods for simulations of Bose gases in a harmonic potential. We usually save 
one thousand or more field configurations at equally spaced times throughout the evolution, which are sufficiently 
separated in time that there is little correlation between samples. 

The temperature analysis is performed after the simulation is completed. The field configurations are stored 
as vectors of basis coefficients, and we make use of efficient computational routines that avoid numerical aliasing 
to transform to and from a real space representation. For each saved field configuration we calculate the vector 
components P{ and Qi, the derivatives in Eqs. ([28l - [3"Tj) , the auxiliary fields given by Eqs. (|40I41|) , and then their 
overlaps with the various functions of the field ?/>(x) as in Eqs. (|54ll58|) . This allows us to calculate the vector 
quantities defined in Eqs. (f70Tl8Tjl that appear in Eq. (|20|). We then form the appropriate dot products of these 
quantities to give a single sample of Dq ■ Xy and T>p ■ Xy. 

These values are calculated for every saved field configuration. These can then be plotted versus time, and any 
initial transients before the field has settled into equilibrium can be identified (typically this a very small fraction of 
the total simulation time.) The quantities T>q ■ X^ and T>p ■ Xy arc numerically distinct from one another before 
equilibrium is reached, however they can be quite noisy and so sometimes it is difficult to tell when they are in 
agreement. We exclude the transients from the final averaging to determine the temperature — typically we discard 
at least the first quarter of the time evolution to ensure accuracy. 

We now illustrate with a set of typical numerical results from applying the method described in this paper to 
the PGPE for trapped Bose gas systems. We solve the PGPE (fT0|) for a harmonic trapping potential V(x) = 
(x 2 + y 2 + 8z 2 )/4 and no additional potential T4d(x). The unit of length is L = (h/mtOx) 1 ^ 2 and energy el = hu x . 
We have a dimensionless energy cutoff of E cut = 31 such that there are 1739 classical modes, and the data shown is 
for C n i = 2000. We begin individual simulations with a randomised initial condition with fixed energy, and evolve 
in dimensionless time for until r = 1200, which is slightly more than 190 radial trap periods. We find that these 
simulations equilibrate very quickly, and use 1000 field samples over the last two-thirds of the evolution for the ergodic 
averaging. The results are shown in Fig. [TJ some of which have previously been reported in Ref. [l5| . 

Figure [TJ a) shows the condensate fraction versus energy for the PGPE system (we remind the reader that energy 
given by Eq. (|12p is a conserved quantity in these calculations) . The condensate fraction is determined by the Penrose- 
Onsager criterion: it is the largest eigenvalue found by diagonalizing the single-particle density matrix formed by 
ergodic averaging as described in [l5| . It looks somewhat different to the result for the full three-dimensional system 
due to the basis cutoff that is present in these simulations. Figure [TJb) shows the distribution of the function T> ■ Xj- 
for two different simulations, using the formulae for both T>q and Dp. The temperature is found by the inverse of the 
mean value of these distributions, and is shown for all simulations in Fig. [TJc). We can see that the results of both 
calculations agree very well, and so we can be confident that firstly the simulations have reached equilibrium, and 
secondly that the numerical implementation of these calculations are free from errors. Figure [TJd) shows the similar 
calculation for the chemical potential. 

One interesting point is that the width of the distributions in Fig.QJb) calculated using the operator T>p are slightly 
narrower than those for T>q . While this is of no consequence conceptually, in practice the narrower the distribution 
the fewer samples are required for an estimate of the mean to a required accuracy. Given the large amount of freedom 
in the choice of the operator T>, it seems quite possible that for particular problems that some choices will be better 
than others. We have found situations where one of these distributions is signficantly narrower than the other. This 
has also been pointed out by Rugh (2l| . However, with no a priori way to estimate the width of the distributions, 
finding the most accurate method of determing the temperature is a matter of trial and error. 

V. CONCLUSIONS 

To investigate the effect of critical fluctuations in Bose gases using the PGPE in the microcanonical ensemble it is 
necessary to have accurate methods of determining the thermodynamic temperature and chemical potential. In this 
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FIG. 1: Data extracted from PGPE simulations in equilibrium, (a) Condensate fraction versus energy of simulation, (b) 
Histograms of the quantity T> ■ Xt on the RHS of Eq. (f4]) for the simulations with E = 10 and E = 11. The bars are for the 
operator T>q, and the dots connected by solid lines for Dp. (c) Temperature versus energy. Plusses: Tq, crosses: Tp. (d) 
Chemical potential versus energy. Plusses: hq, crosses: /ip. All quantities are in dimensionless units as defined in the text. 



paper we have explicitly outlined a method of how to do so using the assumption of ergodicity and the dynamical 
evolution of the PGPE. The method could potentially be applied to other nonlinear classical field Hamiltonians. 
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